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Collective phenomena often feature structured interactions 
commonly conceptualized with complex networks Q. In 
adaptive networks, the interaction structure coevolves with the 
dynamics it supports, yielding a feedback loop that is com¬ 
mon in a variety of complex systems Ena. Understanding 
their asymptotic regimes is a major goal of the study of such 
systems, and an essential prerequisite for applications. In the 
particular case of a dynamic equilibrium, each node in the 
adaptive network undergoes a perpetual change in its state 
and number of connections to other nodes (its degree), while 
a comprehensive set of network measures become stationary. 
A prominent example is the degree distribution, the probabil¬ 
ity distribution of node degrees. For a wide class of adaptive 
networks in dynamic equilibrium, the shapes of stationary de¬ 
gree distributions appear to be insensitive to initial conditions 
in state and topology 110 - not only when taken over the 
whole network (network degree distributions), but also when 
describing ensembles consisting only of nodes of same state 
(ensemble degree distributions). 

While much work on adaptive networks assumes random 
connectivity in the form of Poissonian degree distributions 
El ID, coevolutionary dynamics can generate highly struc¬ 
tured steady-state topologies Em. Analytic expressions for 
the ensuing degree distributions have been so far lacking, and 
their investigation has relied on numerical procedures BE. 
As a consequence, the distributions’ dependency on system 
parameters is difficult to infer and small parameter regions 
with counterintuitive topologies prone to be overlooked. 

Here, we revisit the adaptive contact process in dynamic 
equilibrium ifTTIl . Using a compartmental approach 112], we 
obtain closed-form ensemble degree distributions dependent 
on a single external parameter, and show that a coarse-grained 
understanding of the distributions’ shapes can be obtained 
self-containedly. In particular, the emergence of symmet¬ 
ric ensemble statistics from asymmetric dynamics can be ex¬ 
plained. The framework’s applicability to static networks as 
well as to other coevolutionary dynamics is also discussed. 

Model. — The contact process on an adaptive network mod¬ 
els the spreading of a disease in a population without immu¬ 
nity, but with disease awareness CD. The disease is trans¬ 
mitted along active links that connect infected I-nodes with 


susceptible S-nodes, letting the susceptible end switch to the I- 
state with rate p. Moreover, I-nodes recover to the S-state with 
rate r. Additionally, S-nodes evade infection by retracting ac¬ 
tive links with rate w and rewiring them to randomly selected 
S-nodes. The latter process ties the network’s topological evo¬ 
lution to its state dynamics, yielding an undirected adaptive 
network with constant mean network degree (k). Initial con¬ 
ditions with (k) > 2 should be taken to ensure minimum net¬ 
work connectivity. Ensuing dynamics can be described with 
a low-dimensional pair-approximation (PA) ansatz that tracks 
state correlations among next neighbors CD- 

In the model’s simple active phase, given asw < {k)p—r in 
the PA iflOl . the system reaches independently of initial con¬ 
ditions a dynamic equilibrium characterized by stationary en¬ 
semble degree distributions Pj^{k), A G {S, /}, as well as sta¬ 
tionary values of the fraction [A] of A-nodes and the per-capita 
number [SI] of active links. These values and the form of the 
steady-state Pg i(fc), in particular their first (second) moments 
(A:a) ({k\)), depend only on model parameters. The simple 
active phase is the dominant active regime of the model DD, 
serving as a testbed for our approach laid out in the following. 

Let PA{k) = [A]Pa(A:) be the fraction of nodes of state 
A and degree k in an inhnitely large network with hnite (fc), 
and for a node in state A, denote by /a = [Apfei) average 
fraction of neighbors in the respective other state. Then, the 
coupled state and degree evolution determined by this average 
is given by the master equations 
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In both equations, the hrst, second and third term on the right- 
hand side describe infection, recovery and rewiring, respec¬ 
tively. The third term in the second equation captures the de¬ 
gree gain of S-nodes that are being rewired to. 

Closed-form expressions. — We are interested in the steady 
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state of Eqs. [T] so that all introduced measures are assumed 
stationary for all following considerations. Moreover, the bal¬ 
ance equation p[SI] = r[I] for [S'] and [/] must hold, so that 
with a = w/p, Eqs.[^yield the coupled recurrence relations 

° ^ ^ ^ ^ 

0 = + « ['Ps(fc - 1) - Ps{k)] . (2) 

If a = 0, i.e., in static networks, Eqs. are not independent. 
A fixed network degree distribution 

P{k) = il-[I])Ps{k) + [I]Pi{k) (3) 


degree yields a third independent constraint. With it, the re¬ 
covery rate r enters Eqs.implicitly through the external pa¬ 
rameter [/] as in the static case. But as shown in the following, 
one does not need to undertake the complete solution of Eqs.[^ 
in order to i) infer general properties of Ps,i{k) ii) uncover 
a particular ensemble symmetry iii) considerably reduce the 
search space for self-consistent (fcs,i). 

First moments and symmetry .— Eirstly, we remark that due 
to normalization, Ps{k) and Pi{k) are either identical or in¬ 
tersect at least once. Setting Ps{k) = Pi{k) in Eq. [^reveals 
that there can be at most one such intersection for any choice 
of {ks,i)- Secondly, we see from Eqs. gthat Pi( 0 ) ^ Ps(0) 
iff a ^ 1. Since the distribution dominating the low-degree 
range before the sole intersection possesses the lower mean. 


then determines Psj{k) via the steady-state fraction [J] of I- 
nodes. This fraction must be provided externally, e.g., by the 
model’s PA or simulations of the full system, both observed 
to deliver almost identical values in the simple active phase. 
Note that Eq. holds in general, relating the distributions’ 
moments accordingly when a > 0 and P{k) is not fixed. 

Eor the remainder of this work, the general case a > 0 is 
considered, for which Eqs.j^are in contrast solved by 
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(A;s) ^ (fci)iffa^ 1. (7) 

Hence for a < 1, infection outweighs the rewiring bias to¬ 
wards S-nodes, yielding a higher connectivity of I-nodes cni. 

Erom Eq.|^ follows (fcs) = (^i) iff a = 1- Setting a = 1 
and (fcs) = (fci) in Eqs.|^ we see that then the Pg i(fc) i) coin¬ 
cide ii) are Poissonian with Psj{k) = P{k) = {k)^/k\ 

as in Erdos-Renyi (ER) graphs with same {k) iii) are indepen¬ 
dent of r. Moreover, this is the only choice of a where any of 
assertions i)-iii) hold (see Eqs.0. As laid out in the follow¬ 
ing, a coarse-grained understanding of the Pg i(fc) for a 7 ^ 1 
and beyond Eq.j^can also be obtained algebraically through 
considering their variances and monotonicity. 

Variances .— With the variance tTg j = (fcg j) — (fcgj)^ of 
Pg i(fc) as well as Eqs.|^an d|^ we can relate the moments as 


which in particular implies 

Pg(fc) = Pi(fc + l) 
Pg(fc) = Piik) 
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Eor coevolution with a —)■ 0, Pi(fc) = kPs{k)/{ks) 

andPg(fc) = /o-^( 2 v/(fcs)(A:i))((fcs)(fci))V(fc!)" (Eqs. Q, 
where Io{x) is a modified Bessel function of the first kind. 
This supports the existence, previously conjectured in ||71, of 
a discontinuous transition in the full model from static to coe¬ 
volving steady-state topologies as rewiring is switched on. 

Constraints .— The functional form of the steady-state 
Psi{k) is given by Eqs. whose free parameters 
Pi(0) and (fcgj) can be determined through normaliza¬ 
tion and self-consistency constraints on Ps,i(fc). Obviously 
J2T=o PsAk) = 1 and kPs,iik) = {ks,i) should hold, 

but these constraints are not all independent; from Eq. 1^ 
J2T=okPiik) = (fci) ifEfco-fs(fc) = 1, and from Eq.|^ 
J2T=okPs{k) = (fcs) if Z)fcoEs.i(fc) = 1- Hence normal¬ 
ization implies self-consistency of the first moments, and is 
assumed to be given for all considerations below. 

However, the two normalization constraints obviously do 
not suffice to determine the three free parameters. As (k) = 
(1 — [/])(fcs) -I- [/](fci), imposing a constant mean network 


(te) - (h) = (tly - 1 


( 8 ) 

(9) 


To assess Cg, we set Pi(fc) = fcPg(fc)/(fcg) and P 2 (fc) = 
Ps(fc - 1) = E’i(fc) °fc^f|tfc (Eq. 61 , with Pi, 2 (fc) clearly nor¬ 
malized. Analogously to arguments leading to Eq. yet con¬ 
sidering that Pi(0) = P 2 ( 0 ), we conclude (fci) ^ (fc 2 ) iff 
P 2 (l) ^ Pi(l) andi) (fci) ^ (fc 2 ) iff a(fcg) + 1 ^ (fci) + a. 

Assuming o < 1, it follows that (fcj) > (fcg) (Eq. and 
(fci) > 1 (considering (fc) > 2), so that 1 — a < (l-a)(fci)-f 
a((fci) — (fcs))- Similarly, setting a > 1 yields a — 1 < (a — 
l)(fcs) + (fcs) — (fci)- Hence a(fcg'i -i- 1 ^ 

Inserting Ps(fc) into ii) delivers erg ^ 
moreover ^ (fci) iff o ^ 1 through Eqs.andwe obtain 


(fcg) — (fci). Hence a(fcg) -f 1 ^ (fci) + a iff a ^ 1 
and, with i), ii) X)^o^E’i(fc) ^ J 2 T=okP 2 {k) iff a ^ 1 . 

^ {ks) iff a ^ 1. As 


f^S.i ^ (fcs,i) iff a ^ 1 . 


( 10 ) 


In case of the variance A of P(fc), one concludes from 
Eq.[^that A > (fc) if a > 1. Eurthermore, from Eqs.|^and 
[To]follows |a— 1 | > I (fcg) — (fci) I for a 7 ^ 1 , which with Eqs.[^ 
andj^yields A < (fc) if a < 1 , so that 

A ^ (fc) iff a ^ 1 . 


( 11 ) 
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Thus for a < 1, rewiring actually decreases degree variabil¬ 
ity with respect to ER graphs in both node ensembles (as ob¬ 
served in simulations in jT)) as well as in the overall network. 

Bounds and monotonicity. — From Eqs. □ 0 and 
obtains for a ^ 1 


one 


min(0, a — 1) < {kg) — (ki) < max(0, a — 1) (12) 

with the PA predicting (fcs) — (ki) = a— 1 iflOl . Given a fixed 
(k), EqJ~ 
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; already restricts the range of self-consistent (kg^) 
in Eqs. 1^ independently of the external parameter [/] S [0,1] 
[shaded region in Fig.[2a)], considerably speeding up the nu¬ 
merical solution to the constraint problem. These bounds can 
moreover facilitate heuristic guesses for self-consistent {ksj) 
to make the framework fully self-sufficient. 

For the dispersion indices, 


0’s I 

min(l,a) < ,, ’ , <max(l,a) 
(«s.i) 


(13) 


. 0-^ I 1, if o < 1 

min(l, a) < tw < \ 

{k) I a-I-(a - l)^/4/(fc), if o > 1, 


(14) 


where Eq. [T3|is obtained from Eqs. [8] |9land[T2| and Eq. 
follows from Eqs. [U 12 and TT 


Equations [12 13 and 


set tight bounds for emerging stationary network measures, 
particularly for small a. 




Figure 1: (Color online) (a) Determining self-consistent {fcs,i) for 
a = 2. Solid lines bound the region given by Eq. [T^ the dashed 
line is the actual solution to normalization constraints on Eqs.j^ Im¬ 
posing (k) — 10 further restricts allowed (A:s,i) to the shaded area; 
feeding in [J] = 0.89 as given by the PA forp = r = 1 yields final 
constraint (dash-dotted line). Solution to full constraint problem (red 
circle) is compared to PA prediction (blue square) and simulations 
(green triangle), (b)-(c) Regions (P) in {k),a plane where Ps{k) 
[(b)] and Pi{k) [(c)] peak away from k = 0 for all choices of r, 
compared to boundary of the largest possible (i.e., for r = 0) simple 
active phase in PA (red dashed line). 


From Eqs.|^ it is clear that Ps{k) has one maximum and is 
monotonically decreasing for (fci) < (fcs)~^-Similarly, Pi(A:) 
is monotonically decreasing for (fci) < a and peaks away 
from k = 0 otherwise. For large parameter regions, mono¬ 
tonicity can be assessed by considering how these inequalities 
- together with Eq.[^and fixed (fc) - constrain {kgj). 

It is easy to check that for (fc) > 2, Ps{k) peaks away from 
zero if a < (k) — l/(fc) -F 1 [i.e., in the entire simple active 


phase; Fig. [^b)], whereas Pi(fc) does so for a <m + i)/2 
[Fig. [TJc)]. Numerical investigation of the remaining regions 
yields distributions peaking away from k = 0 throughout the 
simple active phase. 

Comparison to simulations. — With given model parame¬ 
ters (w,p, r, {k)), the external parameter [/] is extracted from 
the PA and self-consistent Psj{k) are generated. As static net¬ 
works, Barabasi-Albert graphs with P{k) ^ k~^ are chosen 
d, and Eqs. 0 and 0 are used for 10^ realizations of P{k) 
[Fig. 0a)]. In the coevolutionary case, Eqs. 0 [Fig. 0d)] or 
one of their limiting cases are employed [a —>^ 0 in Fig. 0b) 
and a = 1 in Fig. 0c)]. To simulate the full system, the 
Gillespie method lfl4i is implemented for network sizes of 
10'* nodes and a runtime of t = 10^, averaging over 10^ re¬ 
alizations. Initial networks are randomly primed with 9 • 10^ 
I-nodes and, in case of subsequent coevolutionary dynamics, 
chosen to be of ER type. 
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Figure 2: (Color online) Stationary statistics of S-nodes (blue cir¬ 
cles/solid lines) and I-nodes (red squares/dashed lines) for (k) = 10, 
r = p = 1 and [/] = 0.89. Main plots: Ps,i{k) in simulations 
(symbols) and model (lines), compared to Poissonian distributions 
of same (fc) (black dotted line). Plot insets: fs,i in simulations; k- 
dependency (symbols) and mean field approximation (lines). Text 
insets: (fcs.i) (<r| i) in simulations, (a) Static Barabasi-Albert graph, 
(b) a = 0.01. (c) a = 1; additionally P{k) for r — 0.2 (black 
crosses, [/] = 0.98) and r — 5 (black pluses, [7] = 0.44). All 
distributions are Poissonian and coincide. (d) a = 3. 

Our framework then delivers excellent predictions if coevo¬ 
lution does not occur on a much faster timescale than disease 
dynamics [Figs.0a)-(c)]. As a is further increased, generated 
distributions start deviating from those obtained from simula¬ 
tions of the full system [Fig. 0d)]. This is because due to 
strong state and degree correlations among next neighbors, 
the mean fields used in Eq. 0 should be replaced by degree- 
dependent expressions fs,i{k), as observed in ||4l for similar 
dynamics (insets of Fig. 0. 

High-degree I-nodes tend to stem from recently infected S- 
nodes, which in turn had accumulated disproportionally many 
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susceptible neigbors due to rapid rewiring. Hence /i(fc) > fi 
and, as the process is cyclic, also fs{k) < /s for large degrees 
k and large a. Conversely, low-degree I-nodes predominantly 
are ’’old”, having had the majority of their susceptible neigh¬ 
bors, but none of their infected neighbors, rewired away. Thus 
fi{k) < fi and fs{k) > fs for small k and large a [plot inset 
of Fig.j^d)]. These degree dependencies let Eqs.[T]underesti- 
mate (overestimate) {ks) ((fci)) for rapid rewiring, accounting 
for deviations observed in Fig.j^d). 

Degree heterogeneity in static Barabasi-Albert graphs does 
not challenge the validity of the mean field approximation for 
/s,i(fc) [plot inset of Fig. a)], corroborating that the mean- 
field breakdown is due to degree correlations induced by rapid 
rewiring im. Note furthermore that in all simulations with 
network coevolution, Eqs. [^andj^as well as the inequalities 
of Eqs. 12 13 and are fulfilled with reasonable accuracy 
[text insets of Figs. ]2m)-(d)]. Simulations moreover confirm 
coinciding Poissonian Psj{k) at a = 1 and for a variety of r, 
i.e., regardless of the steady-state abundancy of the two node 
types [Fig. [^c)]. Remarkably, also other topological mea¬ 
sures of the steady-state network are similar to those in ER 
graphs of same mean degree, so that for a = 1, the asym¬ 
metric coevolutionary dynamics appear to randomize network 
topology. 

Generality .— Poissonian ensemble degree distributions 
also arise in the symmetric coevolutionary voter model m 
with link update. In this model, dynamics are entirely driven 
by active links connecting holders of opinion S with those of 
opinion /; the S-end i) adopts opinion I with rate p or ii) 
rewires the I-end with rate w to randomly selected S-nodes. 
The I-nodes undergo the same dynamics, so that in contrast 
to the adaptive contact process, the model is fully symmetric. 
Then, analogously to the adaptive contact process. 


0 =[5'/]fc 


Pi{k) Psjk) 
ih) {ks) 


a[SI][Ps{k-l)-Ps{k)] 


[57] 


+ ay-y [{k + l)Ps{k -b 1) - kPs{k)] 
(«s) 


(15) 


in steady state, where Eq.[^also holds with swapped indices 
due to the symmetry. Considering the active phase before 
fragmentation transition |l6l[8l, [SI] > 0 must hold. We more¬ 
over assume that the two end nodes of active links are statisti¬ 
cally equivalent, setting Ps{k) = Pi{k) = P{k). Then Eq.fTs] 
becomes (fc -b l)P{k -b 1) = ((fc) -b k)P{k) — {k)P{k — 1) 
and is solved by P{k) = {k)^/k\. Simulations indeed 

reveal coinciding Poissonian Ps,i{k) for large regions of the 
active phase (not shown). 

The proposed framework can be readily applied to other 
two-state coevolutionary models featuring node and link pro¬ 
cesses with constant rates. For dynamics with a larger state 
space, obtaining degree distributions clearly is more involved 
due to the increased number of coupled recurrence relations. 
However, the distributions’ first two moments already provide 
an insightful description of steady-state network topology and 
can often be obtained without the degree distributions at hand. 


Conclusions .— For the adaptive contact process in dy¬ 
namic equilibrium, we use a mean field approximation to ob¬ 
tain closed-form ensemble degree distributions. These are 
parametrized by their first moments which are numerically 
determined through the input of a simple external parameter. 
For small and moderate topological coevolution, ensuing dis¬ 
tributions match very well those observed in the full system, 
while deviations for rapid topology change are explained on 
the basis of a mean-field breakdown. The fraction of rewiring 
and infection rate is identified as the crucial model parame¬ 
ter, allowing for a characterization of the distributions’ shapes 
even without relying on external input. When this fraction 
is smaller (larger) than one, link rewiring is shown to yield 
i) a smaller (larger) mean degree of S-nodes than of I-nodes 
despite its bias ii) less (more) degree heterogeneity than in 
respective ER graphs. When this fraction equals one, we 
show that the asymmetric dynamics is characterized by co¬ 
inciding Poissonian ensemble degree distributions, regardless 
of the value of the recovery rate (see ca for a similar ex¬ 
ample). Apart from explaining these counterintuitive results, 
the method easily generalizes, enabling a quick assessment of 
possible steady-state topologies in adaptive networks. 

Future work could extend beyond second moments the de¬ 
scription of the ensemble degree distributions, and improve 
the simple mean-field assumption used here through account¬ 
ing for state heterogeneity among nodes’ neighbors, in the 
spirit of ||2l . Finally, it has been shown ||5] [TOll that ensem¬ 
ble degree distributions are linked to stationary distributions 
that describe other features of the steady state. Obtaining the 
latter would be another contribution to the study of dynamic 
equilibria in coevolutionary dynamics. 
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